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Abstract 

We explore the gravitational dynamics of falling through planetary interiors. Two trajectory classes are considered: 
a straight cord between two surface points, and the brachistochrone path that minimizes the falling time between two 
points. The times taken to fall along these paths, and the shapes of the brachistochrone paths, are examined for the 
Moon, Mars, Earth, Saturn, and the Sun, based on models of their interiors. A toy model of the internal structure, 
a power-law gravitational field, characterizes the dynamics with one parameter, the exponent of the power-law, with 
values from -2 for a point-mass to -1-1 for a uniform sphere. Smaller celestial bodies behave like a uniform sphere, 
while larger bodies begin to approximate point-masses, consistent with an effective exponent describing their interior 
gravity. 
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This letter considers the dynamics of objects falling along paths constrained to the interiors of planetary bodies, 
where the mass distributions arising from hydrostatic equilibrium lead to non-trivial gravitational fields. Two classes 
of paths are considered: a Euclidean cord connecting two points through the interior, and the brachistochrone path 
that minimizes the time taken to fall between two points. The internal dynamics are calculated for the Moon, Mars, 
Earth, Saturn, and the Sun, and compared to a model in which the internal gravitational field scales as a power-law. 
Paths are assumed to be frictionless and planetary rotation and asphericity is ignored. 

The best-known problem of this type is the gravity tunnel through the Earth [T]. Under the assumption of 
uniform density, the gravitational potential becomes harmonic and the dynamics are those of a simple harmonic 
oscillator. A surprising result of this system is that a cord path takes the same amount of time to fall through, 
42 minutes, independent of the length of the cord. The brachistochrone curve through a uniformly dense sphere is 
that of a hypocycloid, the shape found by rolling a smaller circle inside a larger one [5]. Recently, the gravitational 
dynamics of a falling object inside the Earth were calculated using a realistic internal profile [3], and found to be 
similar to the dynamics under a constant radial field: it is more accurate to assume that the Earth’s gravity is the 
same strength throughout its interior than to assume that the density is constant. This was explained by the fact that 
the internal mass profile of the Earth increases roughly with the 1.97 power of radius, nearly cancelling Newtonian 
gravity. Inspired by that work, this letter extends that analysis to other solar system bodies and attempts to classify 
the dynamics through a single parameter: the scaling exponent of the internal gravitational field. 

We consider the dynamics of a test particle falling through a spherically symmetric gravitational field, where the 
gravitational field scales with radius r as a power-law with an exponent a to some surface field go at surface radius R. 
The mass m enclosed within radius r increases towards M at R with the exponent a -f 2, related to the gravitational 
field through Newton’s law, 


( rp \ ct / r \ Q;+2 

r) rn{r) = M(^-j . ( 1 ) 

There is a certain range of a that is physically meaningful. The value a = 1 corresponds to a uniformly dense 
sphere, whereas for a > 1, density increases with radius and hydrostatic equilibrium is violated. Large values of 
a may describe a hollow shell of mass or a Dyson sphere. The value a = —2 corresponds to Newtonian gravity 
outside a point-mass, thus a < — 2 is unphysical (such fields can be found in region of negative-mass material, or in 
a co-orbiting reference frame around a strong gravitational quadrupole or electric dipole). The density is singular at 
the origin for all a < 1, but avoiding this singularity leads to tension with the constraint of hydrostatic equilibrium. 
The physical values of a are therefore — 2 < a < 1. The special case of a = — 1 corresponds to a logarithmic potential 
and must be treated uniquely, but in practice Q! = —l-|-eisa well-behaved approximation. 

Falling from rest at the surface, the velocity u of a test mass at any radius can be related to that radius by 
conservation of energy. The gravitational potential is found by integrating the field from the zero-point Vo to that 
radius. The zero-point is at the origin for a > — 1 and at infinity otherwise, although ultimately only the difference 
in potential is relevant. 
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The time taken to fall along the diameter from r = i? to r = 0 can be found by integrating the reciprocal velocity 
by the line element, and the time taken to fall to the opposite side is twice this: 
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This closed-form expression is only defined for a > —1 and the value is known from orbital mechanics for a = —2 
[ 4 ], and numerical integration can be used to find the time for the entire range. It is useful to consider this time 
relative to the uniform density (a = 1) value, which is the same as the time taken to make half an orbit along the 

surface: Tr = Tdiameter/7rOver the physical range of a, the relative time Tr is almost perfectly described by 
the linear function = 0.9 -I- 0.1a. 

We now consider a straight cord connecting two surface points, subtending an angle 20 (Figure 1). Making an 
angle 9 with respect to the midpoint between the two points, the radial coordinate of any position along the cord is 
r{6) = i?cos(0) sec{9) and the line element along this path is ds = i?cos(0) sec^{6)d6. The time taken to fall from 
9 = Q to 9 = —0 along a cord can again be found by integrating the reciprocal velocity by the line element. 
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This has a known solution for a = 1 where the cord time is independent of 0. A solution for a = 0 has been 
derived [S]- No general solution for arbitrary a was found, thus numerical analysis is required. The cord time can 
either be found by numerical integration of Equation 5 or the kinematic equation of motion; the latter yields smoother 
results. 

While the cord path minimizes the distance travelled by the falling object, the brachistochrone path minimizes 
the time taken to fall between two points. To find the brachistochrone, the time taken to fall along an arbitrary path 
from point A to point B is written in radial coordinates, following the derivation of the generalized brachistochrone 
[3]. The velocity, known from Equation |3] is left un-expanded for compactness. 
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Because the integrand dt depends only on r and it can be minimized using the Beltrami identity. An additional 
constraint is added, such that the path is horizontal at its deepest point: ^ = 0 when r = Rd, 

C is a constant. This yields the minimization condition for the brachistochrone curve, at which point the power-law 
velocity can be substituted. 
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This differential equation can solved for a relationship between r and 9, and can be integrated to find Tbrac- 
Again, a general solution is lacking but specific cases have been discussed for a = l [2],a = 0 [3],a = —1 [5] , and 
a = —2 [3]. A numerical method [3] is used to calculate the brachistochrone path and time for arbitrary a. 
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Figure 1: Schematic of the coordinate system used in these derivations. A cord path (red) and brachistochrone 
(green) subtend an angle 20 and the brachistochrone reaches a minimum deepest radius at Rd- 

Having derived expressions for falling through a power-law gravitational field along a cord and the brachistochrone, 
we can examine their behavior as a function of a in Figure 2. First considering the time taken to fall through the 
diameter, over the physical range of a, the relative time ranges from Tr = l/%/2 at a = —2 to = 1 at a = 1, with a 
near-linear dependence, and continues to increase in the non-hydrostatic regime. The time taken to fall along a cord 
path is only independent of the surface angle when a = 1. In other cases it interpolates between this orbital half¬ 
period and ^diameter' non-hydrostatic regime, it takes more time to fall along longer paths, as gravitational 

fields weaken near the origin. Brachistochrone curves, which balance the greatest speeds with the shortest paths, 
generally become more vertical as a decreases and there is a greater benefit to a longer, deeper path. For larger a, 
the brachistochrone curves begin to approximate a cord, as the shortest path becomes favored. 

Having considered an ideal case of a power-law gravitational field, we now consider the gravitational dynamics 
inside Solar System bodies. The interiors of the Moon and Mars are given by a piecewise stratified model taken 
from Bills and Rubencam [6]. The interior of the Earth is taken from the Preliminary Reference Earth Model, based 
on seismic data [7]. The interior of Saturn is based on a sixth-order polynomial fit of its internal density profile 
to its gravitational multipole moments, taken from Schubert and Anderson [8], and the interior of the sun is based 
on polytropic hydrostatic equilibrium from Christensen-Dalsgaard [9]. The radial mass and gravity profiles of these 
bodies are seen in Figure 3. 

The relative radial mass profiles of the planetary bodies increase from zero at the origin to unity at the surface, 
although the interpolation is mass-dependent: smaller objects tend to increase more slowly, comparable to a uniformly 
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Figure 2: Left: The time taken to fall to the center of a planetary body as a function of the gravity scaling 
exponent a, normalized by the orbital half-period. The black curve is Equation[4]and the red curve is a kinematic 
numerical solution. Middle. Cord fall times as a function of surface angular separation, normalized to the orbital 
half-period, for several values of a. Right. Brachistochrone paths from an arbitrary surface point, for various 
values of a. 


dense cubic for the moon, while for the larger objects the mass profile is akin to logistic growth as the outer regions 
contribute little mass. Planetary bodies that are in hydrostatic equilibrium typically have a gravitational field that 
is highest in the planetary interior. The Earth, for example, has its strongest gravitational field of 10.8 m/s^ at 
roughly 3000 km below the surface [7]- Lighter objects such as the Moon and Mars have a gradual near-monotone 
increase in gravitational field with radius, consistent with small hydrostatic deviations from uniform density. Larger 
objects, such as the Sun and Saturn, have their mass strongly concentrated towards the center, and have a strong 
sub-surface gravitational peak and a rapid falloff towards zero at the origin. 



Figure 3: Left: The radial mass profiles of the Moon, Mars, Earth, Saturn, and the Sun. For comparison, the 
cases of uniform density (cubic radial dependence) and constant gravity (quadratic radial dependence) are shown. 
Right: The internal gravitational fields of the Moon, Mars, Earth, Saturn, and the Sun. Uniform density (linear 
radial dependence) and inverse-square curves are shown for comparison. The curve for the Sun is truncated so 
that the information in the other curves is visible; it reaches a relative field of 8.6 at a fractional radius of 0.17. 


The physical properties and absolute diameter fall times can be found in Table 1. The time taken to fall through 
a cord between two points is shown in Figure |4j left. All curves are bound by the upper limit of uniform density and 
the lower bound of the point-mass. Curiously, in addition to being the same angular size as seen from Earth, the Sun 
and the Moon would take the same time to fall through. Examining the curves in Figure [d] a trend emerges. For 
smaller objects under hydrostatic equilibrium, gravity deforms the density only slightly, such that the behavior is 
close to that predicted for uniform density, only decreasing by four percent with angular separation, in the case of the 
moon. As the object becomes more massive, the behavior approaches that of a point source as more mass is found 
towards the center. The same trend is evident in the brachistochrone times (Figure 4, center): for smaller objects, 
the behavior is similar to the prediction of the uniform-density sphere, while for larger objects it approaches that of 
the point-mass. The brachistochrone paths (Figure 4, right) are shown for an arbitrary surface angle (corresponding 
to roughly 10,000 km on the Earth), and again this trend is seen, with the curves remaining closer to vertical in 
larger objects. 

We can examine in more detail the line fall times predicted based on the exponents taken from fits to the 
radial mass profiles (Figure [5ll. The exponents are found from single-parameter least-squares fits to the function 
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Table 1: The physical parameters of the celestial bodies, the time taken to fall through their diameters, and 
the effective power-law exponent ot estimated from fits to their mass profiles and from their fall time. 



Radius (km) 

Mass (kg) 

go (m/s2) 

Diameter Time (min) 

Mass a 

Time a 

Moon 

1,734 

7.3 ■ 10^^ 

1.6 

52 

0.6 

0.6 

Mars 

3,390 

6.4 • 10^3 

3.7 

47 

0.4 

0.4 

Earth 

6,380 

6.0 ■ 10^4 

9.8 

38 

0 

0 

Saturn 

60,200 

5.7-1023 

10.4 

61 

-0.8 

-0.9 

Sun 

695,000 

2.0 ■ 103o 

274 

52 

-1.5 

-1.7 



Figure 4: Left. The cord fall times through the Moon, Mars, Earth, Saturn, and the Sun, as a function 
of surface angular separation, normalized to the orbital half-period. Middle: The brachistochrone fall times 
through the Moon, Mars, Earth, Saturn, and the Sun, as a function of surface angular separation, normalized 
to the orbital half-period. Right: Brachistochrone curves through the Moon, Mars, Earth, Saturn, and the Sun, 
with the same surface angular separation throughout. In all figures, the behavior for a uniformly dense object 
(dashed line) and a point-mass (dotted line) are shown for comparison. Note that a brachistochrone curve around 
a point-mass does not exist for a surface angle greater than 27r/3 [4]. 


m(r) = M Fits tend to be biased towards the near-surface trends, but most of time spent falling is also 

near the surface, where velocities are slower. The cord time profiles are qualitatively similar between the exponent 
model and the prediction from the density profile, with deviations for intermediate surface angles. The Sun is better 
described by the dynamics around a point-mass than those found from a fit to its mass profile. To further examine 
the validity of the power-law model, we can compare an “effective exponent” based on the diameter fall time to those 
calculated from power-law fits to the internal mass profiles, seen in Table 1. The exponents from the two measures 
are similar in all cases. Overall this shows that despite the non-monotonicity in the internal gravitational fields of 
these bodies, the power-law model is a good first approximation to their internal dynamics. 

In summary, the time taken to fall through straight and minimal paths through the Moon, Mars, Earth, Saturn, 
and the Sun have been calculated based on their known internal structure and compared to the predictions of a 
simple power-law model, with good qualitative agreement. A more complex general model can be developed to 
better characterize the non-monotonicity of internal gravitational fields. The dynamics of rocky material or dust 
through protoplanetary disks may be characterized by the model developed in this letter, and a natural progression 
of this analysis is to compute a model for the internal dynamics through a galaxy, comparing the predictions of a 
dark matter halo to that of a modified gravitational model. The dynamics through a rotating uniformly dense body 
have been extensively covered by Simoson [10], and an exploration of dynamics through rotating non-uniform bodies 
may be of interest. 

The author acknowledges Grace Dupuis, Chen Karako Argaman, and Kiyoshi Masui for helpful 
comments. 
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Figure 5: The line fall times as a function of surface angle for the five discussed bodies, comparing the prediction 
of a fit exponent to that of the density profile. The case for the point mass is shown for comparison. The y-axis 
in each figure is different. 
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